Benchmarks: Multi-class Logloss

This notebook is all about benchmarking some R code used in this package.

Hardware / Software used:

• Intel i7-4600U
• Compilation flags for C/C++: -O2 -Wall $(DEBUGFLAG) -mtune=core2 (R’s defaults)
• Windows Server 2012 R2
• R 3.3.2 + Intel MKL

Libraries

library(data.table)

data.table 1.10.4
  The fastest way to learn (by data.table authors): https://www.datacamp.com/courses/data-analysis-the-data-table-way
  Documentation: ?data.table, example(data.table) and browseVignettes("data.table")
  Release notes, videos and slides: http://r-datatable.com

library(microbenchmark)
library(Rcpp)
library(ggplot2)
library(plotly)

Attaching package: <U+393C><U+3E31>plotly<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:ggplot2<U+393C><U+3E32>:

    last_plot

The following object is masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:

    filter

The following object is masked from <U+393C><U+3E31>package:graphics<U+393C><U+3E32>:

    layout

# Helper function to print data well in tables
print_well <- function(data, digits = 6) {
  
  # To milliseconds
  data <- data / 1000000
  
  # Sprintf helper
  sprintf_helper <- paste0("%.0", digits, "f")
  
  cat("| Min | 25% | 50% | 75% | Max | Mean |  \n| --: | --: | --: | --: | --: | --: |  \n| ", sprintf(sprintf_helper, min(data)), " | ", sprintf(sprintf_helper, quantile(data, probs = 0.25)), " | ", sprintf(sprintf_helper, median(data)), " | ", sprintf(sprintf_helper, quantile(data, probs = 0.75)), " | ", sprintf(sprintf_helper, max(data)), " | ", sprintf(sprintf_helper, mean(data)), " |  \n", sep = "")
  
  return(data)
  
}
# Test case function
# Arguments renamed to avoid recursive clash
test_case <- function(f, preds, labels, eps) {
  cat("Test case: ", paste(do.call(f, list(preds = preds[1:50],
                                           labels = labels[1:5],
                                           eps = 1e-15)), collapse = ", "), "  \n", sep = "")
}
# Fastest Logloss function
cppFunction("double Lpp_logloss1(NumericVector preds, NumericVector labels, double eps) {
  NumericVector clamped = clamp(eps, preds, 1 - eps);
  NumericVector loggy = -log(1 - clamped);
  double logloss = sum(loggy) / loggy.size();
  return logloss;
}")
# Fastest Transformer function
cppFunction("NumericVector Lpp_vec2mat2vec(NumericVector preds, NumericVector labels) {
  int labels_size = labels.size();
  NumericVector selected(labels_size);
  selected = (preds.size() / labels_size) * seq(0, labels_size - 1);
  selected = selected + labels;
  NumericVector to_return(labels_size);
  to_return = preds[selected];
  return to_return;
}")

Benchmarking Clamped Vector to Logloss

For a 10-class vector of 1,000,000 observations:

• Vector A of length=(1000000 * 10)
• Vector B of length=(1000000) with 10 classes

A = [1:1, 1:2, 1:3, 1:4... 1:10, 2:1, 2:2, 2:3..., 1000000:8, 1000000:9, 1000000:10]
B = [3, 5, 9, 1, 4, 8, 6, ...]

Get the following Vector C, D, and E:

C = [1:4, 2:6, 3:10, 4:2, 5:5, 6:9, 7:7, ...]
D = Clamped C by 1e-15
E = Mean of logloss(D, B)

Initialize data

# Generate random data
set.seed(11111)
data <- runif(10000000, 0, 1)
labels <- round(runif(1000000, 0, 9), digits = 0)
# Background truth example (no clamping though)
array(data[1:50], dim = c(10, 5))

           [,1]       [,2]       [,3]       [,4]       [,5]
 [1,] 0.5014483 0.57219649 0.70236924 0.06440384 0.78924978
 [2,] 0.9702328 0.34292525 0.50889166 0.65780657 0.72449914
 [3,] 0.7876004 0.09627503 0.20268701 0.19973930 0.00886554
 [4,] 0.9022259 0.74235690 0.90706612 0.24664551 0.57348710
 [5,] 0.8141778 0.42274539 0.05441064 0.25997440 0.73034459
 [6,] 0.7998922 0.98402494 0.09045308 0.14869691 0.42458612
 [7,] 0.1158690 0.89258437 0.85759800 0.21375685 0.07684022
 [8,] 0.7171363 0.59541565 0.94780036 0.20598170 0.53318628
 [9,] 0.1020639 0.80531234 0.16479666 0.02075780 0.28521238
[10,] 0.6856938 0.42102379 0.42536097 0.63880218 0.86408083

labels[1:5]

[1] 5 1 1 7 8

data[c(1 + labels[1], 11 + labels[2], 21 + labels[3], 31 + labels[4], 41 + labels[5])]

[1] 0.7998922 0.3429253 0.5088917 0.2059817 0.2852124

-log(1 - data[c(1 + labels[1], 11 + labels[2], 21 + labels[3], 31 + labels[4], 41 + labels[5])])

[1] 1.6088991 0.4199575 0.7110905 0.2306488 0.3357698

mean(-log(1 - data[c(1 + labels[1], 11 + labels[2], 21 + labels[3], 31 + labels[4], 41 + labels[5])]))

[1] 0.6612731

# How many digits for benchmarking in milliseconds
my_digits <- 6L
# How many runs for benchmarking?
my_runs <- 1000L

Benchmarks

# ===== BLOCK 1 =====
faster1 <- function(preds, labels, eps = 1e-15) {
  temp_preds <- Lpp_vec2mat2vec(preds, labels)
  temp_log <- Lpp_logloss1(temp_preds, labels, eps)
  return(temp_log)
}
test_case(faster1, preds = data, labels = labels, eps = 1e-15)

Test case: 0.661273148685013

data1 <- print_well(microbenchmark(faster1(data, labels), times = my_runs)$time, digits = my_digits)

Min	25%	50%	75%	Max	Mean
68.776904	74.052150	77.595495	80.546207	158.218773	79.329577

# ===== BLOCK 2 =====
faster2 <- function(preds, labels, eps = 1e-15) {
  x <- pmin(pmax(preds[((0:(length(labels) - 1)) * (length(preds) / length(labels))) + labels + 1], eps), 1 - eps)
  return(-mean(log(1 - x)))
}
test_case(faster2, preds = data, labels = labels, eps = 1e-15)

Test case: 0.661273148685013

data2 <- print_well(microbenchmark(faster2(data, labels), times = my_runs)$time, digits = my_digits)

Min	25%	50%	75%	Max	Mean
76.919802	83.673872	87.017170	90.827465	168.569882	90.820915

# ===== BLOCK 3 =====
faster3 <- function(preds, labels, eps = 1e-15) {
  x <- pmin(pmax(preds[((0:(length(labels) - 1)) * (length(preds) / length(labels))) + labels + 1], eps), 1 - eps)
  return(-sum(log(1 - x)) / length(labels))
}
test_case(faster3, preds = data, labels = labels, eps = 1e-15)

Test case: 0.661273148685013

data3 <- print_well(microbenchmark(faster3(data, labels), times = my_runs)$time, digits = my_digits)

Min	25%	50%	75%	Max	Mean
76.257605	83.326427	85.514111	88.191601	170.386934	90.041468

# ===== BLOCK 4 =====
cppFunction("double faster4(NumericVector preds, NumericVector labels, double eps) {
  int labels_size = labels.size();
  NumericVector selected(labels_size);
  selected = (preds.size() / labels_size) * seq(0, labels_size - 1);
  selected = selected + labels;
  NumericVector to_return(labels_size);
  to_return = preds[selected];
  NumericVector clamped = clamp(eps, to_return, 1 - eps);
  NumericVector loggy = -(log(1 - clamped));
  double logloss = sum(loggy) / labels_size;
  return logloss;
}")
test_case(faster4, preds = data, labels = labels, eps = 1e-15)

Test case: 0.661273148685013

data4 <- print_well(microbenchmark(faster4(data, labels, eps = 1e-15), times = my_runs)$time, digits = my_digits)

Min	25%	50%	75%	Max	Mean
68.158803	70.104625	79.027087	81.212965	155.654754	77.006338

Summary Results

data_time <- data.table(rbindlist(list(data.frame(Time = data1, Bench = "faster1"),
                                       data.frame(Time = data2, Bench = "faster2"),
                                       data.frame(Time = data3, Bench = "faster3"),
                                       data.frame(Time = data4, Bench = "faster4"))))
data_time <- data_time[, t_mean := mean(Time), by = Bench]
data_time <- data_time[, t_median := median(Time), by = Bench]
data_time$Benchs <- data_time$Bench 
levels(data_time$Benchs) <- paste0("faster", 1:4, "= [", sprintf(paste0("%.0", my_digits, "f"), data_time[, list(min(Time)), by = Bench]$V1), ", ", sprintf(paste0("%.0", my_digits, "f"), data_time[, list(max(Time)), by = Bench]$V1), "], mean=", sprintf(paste0("%.0", my_digits, "f"), data_time[, list(mean(Time)), by = Bench]$V1), ", median=", sprintf(paste0("%.0", my_digits, "f"), data_time[, list(median(Time)), by = Bench]$V1))
my_time <- data_time[, list(min(Time), quantile(Time, probs = 0.25), median(Time), quantile(Time, probs = 0.75), max(Time), mean(Time)), by = Bench]
colnames(my_time) <- c("Function", "Min", "25%", "50%", "75%", "Max", "Mean")
my_time <- my_time[order(Mean, decreasing = FALSE), ]
print(my_time, digits = 6)

Plot Results

ggplotly(ggplot(data = data_time, aes(x = Time)) + geom_histogram(aes(y = ..density..), bins = 20, color = "darkblue", fill = "lightblue") + facet_wrap(~ Benchs, ncol = 2) + geom_vline(aes(xintercept = t_mean), color = "blue", linetype = "dashed", size = 2) + geom_vline(aes(xintercept = t_median), color = "red", linetype = "dashed", size = 2) + labs(x = "Time (Milliseconds)", y = "Density") + theme_bw())

ggplotly(ggplot(data = data_time[, .(Time, Bench)], aes(x = Time, y = ..count.., fill = Bench)) + geom_histogram(aes(y = ..density..), bins = 100, position = "fill") + labs(x = "Time (Milliseconds)", y = "Density") + theme_bw())

ggplotly(ggplot(data = data_time[, .(Time, Bench)], aes(x = Time, y = ..count.., fill = Bench)) + geom_density(position = "fill") + labs(x = "Time (Milliseconds)", y = "Density") + theme_bw())

Scaling Benchmarks

Benchmarker <- function(f, size, runs, digits, name) {
  
  data_runs <- list()
  
  for (i in 1:length(size)) {
    set.seed(11111)
    data <- runif(size[i] * 10, 0, 1)
    labels <- round(runif(size[i], 0, 9), digits = 0)
    cat("  \n  \n## ", name, " run: ",format(size[i], big.mark = ",", scientific = FALSE), " samples (", format(runs[i], big.mark = ",", scientific = FALSE), " times)  \n  \n", sep = "")
    test_case(f, preds = data, labels = labels, eps = 1e-15)
    cat("  \n")
    data_runs[[i]] <- print_well(microbenchmark(f(data, labels, eps = 1e-15), times = runs[i])$time, digits = digits)
    data_runs[[i]] <- data.table(Bench = as.factor(paste0("[", i, "] ", format(size[i], big.mark = ",", scientific = FALSE))), Function = as.factor(name), Time = data_runs[[i]])
    gc(verbose = FALSE)
  }
  
  return(data_runs)
  
}
bench_size <- c(100, 1000, 10000, 100000, 1000000, 10000000)
bench_runs <- c(10000, 5000, 1000, 500, 100, 50)

run1 <- Benchmarker(faster3, bench_size, bench_runs, my_digits, "Pure R")

Pure R run: 100 samples (10,000 times)

Test case: 1.35212558439888

Min	25%	50%	75%	Max	Mean
0.022808	0.025469	0.027370	0.043335	4.638802	0.036441

Pure R run: 1,000 samples (5,000 times)

Test case: 1.11833540263107

Min	25%	50%	75%	Max	Mean
0.073366	0.076787	0.083630	0.091613	3.864085	0.088550

Pure R run: 10,000 samples (1,000 times)

Test case: 0.797706755821707

Min	25%	50%	75%	Max	Mean
0.586930	0.595673	0.670181	0.683865	3.591907	0.670329

Pure R run: 100,000 samples (500 times)

Test case: 1.90925095682509

Min	25%	50%	75%	Max	Mean
6.245258	6.728032	7.275997	7.517289	29.690539	7.664803

Pure R run: 1,000,000 samples (100 times)

Test case: 0.661273148685013

Min	25%	50%	75%	Max	Mean
64.671813	67.926919	71.717638	75.255851	151.530276	74.316310

Pure R run: 10,000,000 samples (50 times)

Test case: 0.966169328545996

Min	25%	50%	75%	Max	Mean
741.884853	786.242472	809.841423	821.969478	993.530082	812.830624

run2 <- Benchmarker(faster4, bench_size, bench_runs, my_digits, "Rcpp")

Rcpp run: 100 samples (10,000 times)

Test case: 1.35212558439888

Min	25%	50%	75%	Max	Mean
0.006082	0.006463	0.006843	0.007223	0.068045	0.007184

Rcpp run: 1,000 samples (5,000 times)

Test case: 1.11833540263107

Min	25%	50%	75%	Max	Mean
0.051318	0.055120	0.060062	0.067284	0.210595	0.064723

Rcpp run: 10,000 samples (1,000 times)

Test case: 0.797706755821707

Min	25%	50%	75%	Max	Mean
0.510143	0.515369	0.583890	0.593773	1.034731	0.574722

Rcpp run: 100,000 samples (500 times)

Test case: 1.90925095682509

Min	25%	50%	75%	Max	Mean
5.590283	5.684082	5.804680	6.427818	12.891180	6.154242

Rcpp run: 1,000,000 samples (100 times)

Test case: 0.661273148685013

Min	25%	50%	75%	Max	Mean
57.931618	59.823270	62.286268	66.095328	138.984261	63.797930

Rcpp run: 10,000,000 samples (50 times)

Test case: 0.966169328545996

Min	25%	50%	75%	Max	Mean
673.112510	692.100884	695.991768	706.077447	800.154031	702.179116

Scaling Results

run1_all <- rbindlist(run1)
run2_all <- rbindlist(run2)
run_all <- rbind(run1_all, run2_all)
run_all$Repeats <- rep(inverse.rle(list(lengths = bench_runs, values = bench_size)), 2)
run_all$MilObs <- (1000 / run_all$Time) * run_all$Repeats / 1000000
run_time <- run_all[, list(quantile(Time, probs = 0.05), median(Time), quantile(Time, probs = 0.95), mean(Time)), by = list(Function, Bench)]
colnames(run_time) <- c("Function", "Benchmark", "5%", "50%", "95%", "Mean")
run_time$`Mil.Obs/s` <- (1000 / run_time$Mean) * bench_size / 1000000
run_time$`5%` <- format(run_time$`5%`, digits = 6, scientific = FALSE)
run_time$`50%` <- format(run_time$`50%`, digits = 6, scientific = FALSE)
run_time$`95%` <- format(run_time$`95%`, digits = 6, scientific = FALSE)
run_time$Mean <- format(run_time$Mean, digits = 6, scientific = FALSE)
run_time$`Mil.Obs/s` <- format(run_time$`Mil.Obs/s`, digits = 6, scientific = FALSE)
print(run_time[1:(nrow(run_time) / 2)])

print(run_time[(nrow(run_time) / 2 + 1):nrow(run_time)])

ggplot(run_all, aes(x = Bench, y = Time, color = Function, fill = Bench)) + geom_boxplot() + scale_y_log10(labels = scales::comma, breaks = c(0.01, 0.1, 1, 10, 100, 1000, 10000)) + stat_summary(fun.y = mean, geom = "line", aes(group = Function)) + stat_summary(fun.y = mean, geom = "point", aes(group = Function)) + labs(x = "Benchmark", y = "Time (Milliseconds)") + theme_bw()

ggplot(run_all, aes(x = Bench, y = MilObs, color = Function, fill = Bench)) + geom_boxplot() + scale_y_log10(labels = scales::comma, breaks = c(1, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25), limits = c(1, NA)) + stat_summary(fun.y = mean, geom = "line", aes(group = Function)) + stat_summary(fun.y = mean, geom = "point", aes(group = Function)) + labs(x = "Benchmark", y = "Throughput (Million Obs./s)") + theme_bw()

ggplot(data = run_all, aes(x = MilObs, color = Function, fill = Function, group = Function)) + coord_flip() + stat_ecdf(aes(ymin = ..y.., ymax = 1), alpha = 0.5, geom = "ribbon") + stat_ecdf(geom = "line", size = 2, alpha = 0.75, pad = FALSE) + labs(x = "Throughput (Million Obs./s)", y = "Percentile") + facet_wrap(~ Bench, dir = "h", ncol = 2, scales = "free") + theme_bw()